Optimized Monte Carlo Method for glasses 

L. A. Fernandez 1,3 '* V. Martin-Mayor 1 ' 3 '^ and P. Verrocchio 2,3 ' 4 * 

1 Dep. de Ffsica Teorica I, U. Complutense, 28040 Madrid, Spain. 
2 Dip. di Fisica, U. di Trento, 38050 Povo, Trento, Italy. 
3 Institute de Biocomputacion y Ffsica de Sistemas Complejos (BIFI). Zaragoza, Spain. 
4 Research center Soft INFM-CNR c/o U. di Trento 38050 Povo, Trento, Italy. 

February 6, 2008 



Abstract 

A new Monte Carlo algorithm is introduced for the simulation of supercooled liquids and glass 
formers, and tested in two model glasses. The algorithm is shown to thermalize well below the 
Mode Coupling temperature and to outperform other optimized Monte Carlo methods. Using the 
algorithm, we obtain finite size effects in the specific heat. This effect points to the existence of a 
large correlation length measurable in equal time correlation functions. 
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1. Introduction 

The lack of structural or thermodynamic changes at the glass transition ,1 is a major problem for its 
investigation. The only standard feature, as compared with second-order phase transitions [2], is the 
dramatic dynamical slowing down at the critical temperature. A fairly standard mechanism for slow 
dynamics in an homogeneous system at finite temperature is the divergence of a correlation length 
(critical slowing down 2 ). Slowness arises from the need of configurational changes to propagate over 
increasingly large regions (the critical origin of the Mode Coupling singularity 3 has been recently 
recognized 4 ). It has been recently proposed [5] to study this growing lengthscale in glassformers 
through the finite size behaviour of small systems 2 . Note that experiments in films and nanoporesOEj 
show that the glass transition changes in samples with one or more dimensions of nanometric scale. In 
particular, the specific heat is most sensitive to the size of the confining pore when temperatures are 
close to the glass transition [5]. 

Numerical simulations are an important tool for the study of the glass transition. Their worse 
drawback is the shortness of the times that may be simulated in today computers (roughly speaking, 
microseconds). As a consequence, the computer model goes out of equilibrium by the Mode-Coupling 
temperature, T mc , rather than the actual glass temperature, T c . To approach T c , one may resort to 
optimized Monte Carlo (MC) methods [9"l ll0l l5]. namely methods implementing unphysical dynamical 
rules that strongly reduce the equilibration times. When thermalization is achieved, optimized MC 
allow to study equilibrium mean values and their temperature (or pressure) derivatives, although the 
purely dynamic features of these methods are interesting on their own right 5 . 

Here, we give the first full description of the local swap algorithm 0. We compare the performance 
of the local swap dynamics with the standard MC and with the microcanonical algorithm JO]. We 
conclude that local swap yield equilibrium data at temperatures where the microcanonical algorithm 
no longer thermalizes. Finally, we address the fishy issue of estimating the specific heat in a metastable 
liquid state. We give here details on the strategy followed in where tiny but clearly measurable 
finite-size effects were observed in the specific heat. 

2. Models and observables 

We consider two similar models of fragile glass formers, namely binary mixtures of soft spheres. The 
first model, extensively studied in Ref. |5J , is a 50% mixture of particles interacting through the pair 
potential V a p(r) — e[(<r a + ap)/r\ 12 + C a p, where a,f3 = A,B, with a cutoff at r c = \/3ao- The choice 
a a = 1-2<7_b hampers crystallization, as compared with the oa — cb model. We impose (2(7a) 3 + 
2(tr J 4 + o\b) 3 + (2<tb) 3 = 4<7q where ao is the unit length. Constants C a p are chosen to ensure continuity 
at r c . The simulations are at constant volume, with particle density fixed to c^ 3 and temperatures in 
the range [0.897T mc , 10.792T mc ]. We use periodic boundary conditions on a box of size L in systems 
with N = 512,1024,2048 and 4096 particles. For argon parameters, er = 3.4A, e/k B = 120K and 
T mc = 26.4K. 

In the second model the choice <jb = I-^cta is made. Naming x — r/(a a + ap), the pair 

potential in units e = 1 is V(x) = x~ 12 +x— 13/12 12 / 13 if x < 12 1 / 13 and zero otherwise (thus, the cut-off 
distance depends on the type of interacting species). We study density 1.08cr,^ 3 , as in Refs. [Ill IIP) . 

Since the potential energy per particle, e shows the slowest excitations E] > we shall focus here 
only in this observable, the internal energy being T^k^T + (e) (for other quantities, see Ref. [S]). The 
constant- volume specific heat, C v , is: 

e= 2^E^-^)' C « = I + S [< e2 > " < e > 2 ] (units e = fc B = l). (1) 
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Figure 1: The local swap move. The A particles are depicted by full symbols, while B particles are 
open symbols. The (say) B particle picked randomly (signalled by a central mark) may be swapped 
with any one of the three A particles inside the sphere shown in the left part of the plot. After the 
swap (right), only two B particles could be exchanged with the picked B particle. 



3. The local swap Monte Carlo algorithm 

The Grigera-Parisi swap algorithm 9 consists in picking randomly a pair of particles of distinct type, A 
and .B, try to exchange their positions, and accepting this move with probability min{l, e~ AE / T } (AE 
is the total potential energy change produced by the swap). If combined with standard MC, it is very 
effective in reducing the equilibration time. Nevertheless, there is a caveat: the acceptance of the swap 
move is very small, and it significantly decreases when the number of particles increases. Indeed, the 
closer the swaped particles are, the larger the swap acceptance becomes (this is a huge effect). With 
large systems, it is higly improbable to pick for the swap neighboring particles. To cure this problem, 
we have proposed the local swap. 

In the local swap, the elementary MC step is either (with probability p) a single-particle displacement 
attempt or (with probability 1 — p) an attempt to swap particles. Therefore, for p — 1 the algorithm 
reduces to standard MC. From here on we call local swap to the algorithm with p = 0.5. The time 
unit to is N/p elementary steps. Our swap consists in picking a particle at random (the picked particle) 
and trying to interchange its position with that of a particle of opposite type (the swapped particle) , 
chosen at random among those at distance smaller than 0.6r c (see Fig. Yet, there is an intrinsic 
lack of symmetry. Indeed, let JV id be the number of swappable particles around the picked particle 
in its original position, and N new the number of swappable particles for the picked particle in its final 
position. The probability of choosing the swaped particle is l/7V id in the original configuration, while it 
would be l/./V new in the final configuration. Detailed balance (see e.g. 0) holds only if this asymmetry 
is incorporated in the probability of accepting the swap: 

. M N \d _A£/Ti , n \ 

Paccept swap = nunjl, — — e 1 \. (2) 

1 * new 

The acceptance of the local swap is independent of the number of particles, and larger by a factor 
of 10 than for the original swap algorithm already for N = 1024. For the model cfa/^b — 1-2, 
the acceptance varies from 0.74% at 0.9T mc up to 6% at 2X mc . For the (JaI^b = 1-4 is much smaller, 
actually of the order of 8 x 10~ 6 at T = 0.83T mc , as could be guessed by the disparity in particle 
diameters. 

The performance of the local swap below T mc is far superior to the standard MC (FigEJl and to 
the microcanonical method of Ref. (FigEJl- Furthermore, for both the <ja/&b = 1.2 and the 
ca/c_b = 1.4 models, local swap finds a crystallization phase transition, to highly disordered crystals 
on the bec family, not reported on previous studies [111 ITU| . 
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Figure 2: Comparison of the performance of the local swap (right) and standard Monte Carlo (left), 
as shown by the Monte Carlo history of the potential energy per particle (data points are the average 
of 10 4 io succesive steps), for the <jaI&b = 1-2 model at T = 0.897T mc and 1024 particles. The two 
standard Monte Carlo runs had as starting points thermalizcd configurations obtained with local swap. 
In our time window, the two standard simulations do not explore the same energy range. Instead the 
single local swap simulation explores the full energy range. 
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Figure 3: Specific heat of the (ja/vb — 1-4 model with 128 particles, as a function of temperature, 
as obtained with local swap (10 9 io steps) and with the microcanonical method [10) (error estimates 
were not provided in Ref. [10)). Local swap produces a significantly larger estimate of the specific heat 
at our lowest simulated temperatures, signalling better sampling of configuration space and indicating 
that the microcanonical method is unable to thermalize at such low temperatures. 



4 



2.05 
w 2.00 
1.95 
2.05 
w 2.00 
1.95 
2.05 
a) 2.00 
1.95 





4.5 
4.4 
4.3 
4.2 
4.1 



lxlO 7 2x10' 
t/t 



3x10 



7 





i 


1 

r=o.897r mc , 












i 








m 


r=o.92ir mc x " 






10 


15 

L/v 



Figure 4: (Left) Examples from the 100 generated (local swap) Monte Carlo histories for the potential 
energy density of the gaI^b = 1-2 model with 1024 particles at T — 0.897T mc . One may clearly 
distinguish the metastable liquid from the crystallizing system, that starts as a sharp energy drop in 
all three runs. (Right) Specific heat for the model (JaI g b — 1-2 for sistems with 512, 1024, 2048 and 
4096 particles, at T — 0.921T mc (empty symbols) and T — 0.897T mc (full symbols), versus the size 
of the simulation box (the horizontal line is a fit of the T — 0.921T mc data to a constant value, with 
X 2 /d.o.f. = 4.89/3). At T — 0.897T mc , the specific heat increases with system size, while the energy 
density (not shown), is independent of the number of particles, within our accuracy. 



4. Finite size effects 

Even if not previously known, the liquid state is metastable in our soft-sphere models. Thus, the 
question arises of how to study the thermodynamic properties of a metastable phase. The underlying 
assumption is that the equilibration time for the metastable liquid phase is much smaller than the 
crystallization time. Our strategy has been to run several MC runs (up to 400 at the lowest temper- 
atures .5 ). On each run, the equilibrated metastable liquid is neatly separated from the crystallizing 
system (see Fig. QJ-left). In the analysis we only consider histories whose metastable liquid part was 
selfconsistenly found to be longer than 100 exponential autocorrelation times (to avoid bias, we also 
discarded some 20 autocorrelation times before crystallization). Note that the central MC history in 
Fig- EJ-left does not meet this criterium (at this low temperature the autocorrelation time was 10 5 to)- 
Then, one uses Eq. , with the metastable liquid part of the acceptable histories to obtain the results 
in Figure fright (data from Ref. |S]). 

5. Conclusions 

We have given the first full description of the local swap algorithm py. We compare the performance of 
the local swap dynamics with the standard MC and with the microcanonical algorithm |10) . concluding 
that local swap yield equilibrium data at temperatures where the microcanonical algorithm no longer 
thermalizes. Furthermore, local swaps finds crystal states not reported in previous work |10| . We have 
shown how to estimate the specific heat in a metastable liquid state, finding at low temperatures tiny 
but clearly measurable finite-size effects 
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